# -*- coding: utf-8 -*-
# Code sous licence CeCill
# http://www.cecill.info/licences.en.html
# This is a license compatible with the GNU GPL license
# and that conforms to the principles of French law.
# Authors:
#       Nicolas Bruyant, CNRS,LNCMI
#       Benjamin Edwards, CNRS,LNCMI
"""
Import Filter
=============

Utility functions to import datas and calculate field from pick-up datas

Functions index:
    loadfile
    field_integration
    field_adjust

TODO:
----
    Improve loadfile function to allow column naming
    Add logging support

"""
__version__ = "Version 0.1, SVN-$Revision: 18 $"
"""
Revision allow a good tracking of data processing  errors
Refers to SVN revisions: DO NOT edit manually
"""

import numpy as np

from scipy import integrate,mean,arange,polyfit,poly1d



def loadfileFFC(filename,sample_rate=500E3,nb_col=5):
    """
    Load raw data file from a file of type binary FFC

    Parameters
    ==========
        Filename,
        Sample rate of the datas,
        Number of column of the datafile

    Returns
    =======
        Numpy array containing datas

    Examples
    ========
    >>> inputfile='/Users/nico/Documents/Manips/Superadiance/Documents/Toulouse experiment/Data/jva004'
    >>> temp=loadfile(inputfile,10E6,4)

    """
    with open(str(filename),'rb') as fd:
        temp=np.fromfile(fd,np.float32,-1,"")
    temp.resize(len(temp)/nb_col,nb_col)
    return temp


def loadfileGage(filename,sample_rate=1E6):
    """
    Load raw data file from a file of type binary GageScope

    Parameters
    ==========
        Filename,
        Sample rate of the datas (unused)

    Returns
    =======
        Numpy array containing datas

    Examples
    ========
    >>> filename="/Users/Nicolas/Documents/Projets/TDO/post traitement/tracking freq/s120725d.sig"
    >>> temp=loadfile(inputfile,10E6,4)

    """
    with open(str(filename),'rb') as fd:
        fd.seek(512)   # Skip useless header
        temp=np.fromfile(fd,np.int16,-1,"")
    
    #skip header    
    #temp = temp[range(256,np.size(temp),1)]
    
    return temp


def field_integration(dbdt_data,sample_rate=500E3,method='all',
                      offset_calculation_length=100):
    """
    Integrate pick-up signal for the field

    Parameters
    ==========
        dbdt datas,
        Sample rate,
        Method of offset substraction (Optional):
            -"all": Mean of the whole dataset (default)
            -"start": Mean of the first N samples
            -"end": Mean of the last N samples
        offset_calculation_length: number of sample for offset calculation

    Returns
    ======
        Integrated pick-up signal

    """
    #test for saturation of PU
    absPU=np.asarray(abs(dbdt_data))
    PUmax=np.amax(absPU)
    SaturationRatio=np.count_nonzero(absPU==PUmax)/len(absPU)
    print "saturation ratio:", SaturationRatio
    
    #TODO make it work ...
    if  SaturationRatio > 1e-5:
        print "Saturated pick-UP"
    else:
        print "PU OK"
    if method == 'end':
        offset=np.mean(dbdt_data[-offset_calculation_length:])
    elif method == 'start':
        offset = np.mean(dbdt_data[0:offset_calculation_length])
    else:
        offset = np.mean(dbdt_data)


    temps=np.arange(0,len(dbdt_data))*1/sample_rate
    field=integrate.cumtrapz(dbdt_data-offset,temps)

    return field


def field_adjust(champ, method='end'):
    """
    Adjust the calculated field signal to remove residual offset

    Parameters
    ==========
        champ:  array
            Integrated pick-up signal
                Method of offset substraction:
            -"start": linear fit on the begining of the pulse (default)
            -"end": linear fit on the end of the pulse

    Returns
    =======
        Corrected signal

    TODO
    ====
    utiliser le fit automatique sur toute la plage
    """
    champ=np.asarray(champ)    
    # threshold pour un champ superieur à 1% du max
    threshold = np.max(np.abs(champ))/100

    for i in range(len(champ)):
        if abs(champ[i]) > threshold:
            break
    else:
        print "pb de seuilavant le pulse"

    #on se place 10% ech avant le threshold
    imin = int(i*0.9)
    for i in range(i, len(champ)):
        if abs(champ[i]) < threshold:
            break
    else:
        print "pb de seuil après le pulse"
    #on se place 50% ech après le threshold
    imax = int(i*1.5)

    if method == 'end':
        x2 = arange(imax, len(champ))
        y2 = champ[imax:]
        m = polyfit(x2, y2, 1)
    else:
        x1 = arange(0, imin)
        y1 = champ[:imin]
        m = polyfit(x1, y1, 1)

    offset = poly1d(m)
    champ_corr = champ - offset(arange(0, len(champ)))
    print "max field= %e T" % max(champ_corr)
    return champ_corr
